### ### ### ### ### Replication script cluster analysis and figures ### #### ### ### ###

### load data
load("XXXX/Replication_Data_ClusterAnalysis.RData")

### load library
library(ggplot2)
library(RColorBrewer)
library(gplots)

########################################################################################
################################# Cluster Analysis #####################################

mydata <- na.omit(welfare_indicators) # listwise deletion of missing
rownames(mydata) <- mydata$GEO
mydata$GEO <- NULL
mydata <- scale(mydata) # standardize variables

# Determine number of clusters
wss <- (nrow(mydata)-1)*sum(apply(mydata,2,var))
for (i in 2:15) wss[i] <- sum(kmeans(mydata, 
                                     centers=i)$withinss)
plot(1:15, wss, type="b", xlab="Number of Clusters",
     ylab="Within groups sum of squares")

# Ward Hierarchical Clustering
d <- dist(mydata, method = "euclidean") # distance matrix
fit <- hclust(d, method="ward.D2") 
plot(fit) # display dendogram

# draw dendogram with red borders around the 4 clusters 
rect.hclust(fit, k=4, border="red")

########################################################################################
####################################### Figure 1 #######################################

datascaled <- as.matrix(mydata)
groups <- cutree(fit, k=4) 

cols <- brewer.pal(max(groups), "PuBuGn")

dev.off()
heatmap.2(datascaled, scale='none', 
          col=colorRampPalette(brewer.pal(8, "Blues"))(25),
          cexRow=1.5,cexCol=1.5,margins=c(12,10),trace="none", srtCol=45, 
          RowSideColors=cols[groups])


########################################################################################
####################################### Figure 2 #######################################

dev.off()
jpeg("welfare_map_blues.jpeg", width = 8, height = 10, units = 'in', res = 300)
ggplot(data = welfare_df_map) +
  geom_sf(data= welfare_df_map, mapping = aes(fill = as.factor(groups)), color= "white") +
  scale_fill_brewer(palette="PuBuGn", 
                    name= "Welfare Cluster", 
                    breaks=c(1, 2, 3, 4), 
                    labels= c("Continental welfare cluster", 
                              "Southern-Mixed welfare cluster", 
                              "Eastern European welfare cluster", 
                              "Nordic-Atlantic welfare cluster ")) +
  coord_sf(xlim = c(-21, 40), ylim = c(33, 71), expand = FALSE) +
  theme(panel.grid.major = element_line(color = gray(.9),
                                        linetype = "dashed", size = 0.5),
        panel.background = element_rect(fill = "white"))

############################# Figure 3 was made in Gephi ###############################
